/*==================================================
Project:       Targeting Social Programs
Authors:       Diether W. Beuermann
               Bridget Hoffmann        
               Marco Stampini 
               David L. Vargas
               Diego Vera-Cossio
----------------------------------------------------
Creation Date:    May 2023
Modification Date:   
Do-file version:    01
References:          
Output:             
==================================================*/

/* Regular PMT model using all SISBEN varibles and SISBEN income
 This will serve as baseline for comparison		*/


/*==================================================
            0: Program set up
==================================================*/
*Written on STATA 17
drop _all

** source dir
cd "${dir5r}" // graph dir 

* ssc install winsor2

/*==================================================
            1: load and transformations
==================================================*/

*----------- 1.0.1 Load CRRA welfare FUN
qui do "${dir6r}/ados/bstrap_pmt.ado" // load bstap FUN
qui do "${dir6r}/ados/bstrap_pmt_mh.ado" // load bstap FUN
qui do "${dir6r}/ados/baseline_pmt.ado" // load baseline FUN
qui do "${dir6r}/ados/expanded_pmt.ado" // load expanded FUN

*----------  1.1. Data prep:
use "${dir3r}/01_survey/survey_targeting_r1_long.dta", clear

* new interactions 
cap drop d_loss_illness_inc
cap drop d_reco_ilness_inc
gen d_loss_illness_inc = illness_incidence * d_loss_ilness
gen d_reco_ilness_inc = lag_illness_incidence * d_reco_ilness
lab var d_loss_illness_inc	"Loss: Illness Incidence"
lab var d_reco_ilness_inc	"Recovery: Illness Incidence"


cap drop d_loss_cut_remittance_inc
cap drop d_reco_cut_remittance_inc
gen d_loss_cut_remittance_inc = cut_remittance_incidence * d_loss_cut_remittance
gen d_reco_cut_remittance_inc = lag_cut_remittance_incidence * d_reco_cut_remittance
lab var d_loss_cut_remittance_inc "Loss: Remittances Cut Incidence"
lab var d_reco_cut_remittance_inc "Recovery: Remittances Cut Incidence"

* Globals with shocks
glo loss d_loss_ilness d_loss_cut_remittance death divorce bankrupt nat_shock
glo loss_rec      d_loss_ilness d_reco_ilness d_loss_cut_remittance d_reco_cut_remittance death divorce bankrupt nat_shock 
glo loss_inc      illness_incidence cut_remittance_incidence death divorce bankrupt nat_shock 
glo loss_rec_inc  d_loss_illness_inc d_reco_ilness_inc d_loss_cut_remittance_inc d_reco_cut_remittance_inc death divorce bankrupt nat_shock 


*Globals with updted survey variables
glo hh_char_srv urban age_feb20_srv prop_kids_srv n_members_srv elementary_srv highschool_srv tertiary_edu_srv children_srv living_couple
glo assets_srv own_washing_machine_srv own_tractor_srv own_moto_srv own_car_srv own_PC_srv own_fridge_srv
glo hh_srv  wall_finished_srv floor_finished_srv cookpower_connected_srv wc_flush_srv kitchen_srv electricity_srv running_water_srv trash_service_srv 
glo labour_srv formal_work_srv informal_work_srv

// Including labour status from SISBEN
glo pmt_1_srv  $hh_char_srv $assets_srv $hh_srv $labour_srv

** make sure all labels are OK 
lab var illness_incidence 			"Illness Incidence"
lab var cut_remittance_incidence 	"Remittances Cut Incidence"
lab var death						"Death"
lab var divorce					    "Separation of Spouses"
lab var bankrupt					"Bankruptcy or Business Closure"
lab var nat_shock					"Natural Disaster or Fire"
lab var job_loss                   "Loss: Job lost at any point of the year"
lab var job_loss_formal            "Loss: Job lost formal at any point of the year"
lab var job_loss_informal          "Loss: Job lost informal at any point of the year"

*----------  1.1.2 Estimate PMT for survey data:


//============== Dynamic estimation 

// -------- OLS
local nreps=1000

//define poverty line in terms of Colombian Pesitos
local base=1.90*1515.1
// 2019 PPP factor https://data.worldbank.org/indicator/PA.NUS.PRVT.PP?locations=CO

** new (go back to this in a bit)

// only negative changes in assets
baseline_pmt, pmtvars(l_pp_inc_srv $pmt_1_srv)  predvars($hh_char_srv $hh_srv $labour_srv) negonly($assets_srv) reps(`nreps')  base(`base')
mat amat = (1\1\1), r(estt)
mat estmat0 = amat

// labour updates
baseline_pmt, pmtvars(l_pp_inc_srv $pmt_1_srv)  predvars($hh_char_srv $hh_srv $assets_srv) negonly() reps(`nreps')  base(`base')
mat amat = (2\2\2), r(estt)
mat estmat0 = estmat0 \ amat

// labour updates and assets update 
baseline_pmt, pmtvars(l_pp_inc_srv $pmt_1_srv)  predvars($hh_char_srv $hh_srv) negonly($assets_srv) reps(`nreps')  base(`base')
mat amat = (3\3\3), r(estt)
mat estmat0 = estmat0 \ amat


** send to stata
cap frame create points
frame change points 

** append matrices

** dynamic points
// OLS
drop _all 
svmat estmat0, names(col)
gen type_m = "OLS"
g model_n = c1 
replace type_m =  "negassets"  if model_n == 1
replace type_m =  "labour"   if model_n == 2
replace type_m =  "labour-assets"  if model_n == 3
replace model_n = 12 + c1 if c1 != .
tempfile dfdyn
save `dfdyn', replace


cap drop base _merge 
drop c1

foreach v in covered inclusion_err inclusion_ok exclusion_err exclusion_ok zero epov {
    replace `v' = `v' * 100
    replace `v'_uci=`v'_uci*100
    replace `v'_lci=`v'_lci*100
}

 // PPP values 

 gen hh_benefits_lcu = hh_benefits
 replace hh_benefits = hh_benefits / 1515.1 
 replace hh_benefits_uci= hh_benefits_uci/1515.1
 replace hh_benefits_lci= hh_benefits_lci/1515.1
 // 2019 PPP factor https://data.worldbank.org/indicator/PA.NUS.PRVT.PP?locations=CO

append using "${dir3r}/03_auxiliar/Boostrap_estimates.dta" 

foreach v in covered hh_benefits inclusion_err inclusion_ok exclusion_err exclusion_ok crra crra_l crra_u zero epov covered_uci covered_lci inclusion_err_uci inclusion_err_lci inclusion_ok_uci inclusion_ok_lci exclusion_err_uci exclusion_err_lci exclusion_ok_uci exclusion_ok_lci crra_uci crra_lci crra_l_uci crra_l_lci crra_u_uci crra_u_lci hh_benefits_uci hh_benefits_lci budget_nat_uci budget_nat_lci zero_uci zero_lci epov_uci epov_lci {
    sum `v' if year == 2019 & model_n == 1
    replace `v' = r(mean) if year == 2019 & model_n != 1 
}

// store data for later use 
*save "${dir3r}/03_auxiliar/Boostrap_estimates_lassoRob.dta" , replace
 
/*==================================================
             Make graphs
==================================================*/

foreach m_type in "pmt_dyn"  {

    *** Gain and losses
    #delimit ;
    tw (connected exclusion_err year if type_m ==  "Baseline", 
            mcolor(red%80) msym(S) lcolor(red%80) )
			(rcap exclusion_err_uci exclusion_err_lci year if type_m ==  "Baseline", 
            lcolor(red%80) )
        (connected exclusion_err year if type_m ==  "labour", 
            lcolor( ltblue)  mcolor( ltblue) msym(0) lpattern(-))
			(rcap exclusion_err_uci exclusion_err_lci year if type_m ==  "labour", 
            lcolor( ltblue) lpattern(-))
        (connected exclusion_err year if model_n ==  6 & type_m ==  "OLS" , 
            mcolor(black%80) msym(D) lcolor(black%80) ) 
		(rcap exclusion_err_uci exclusion_err_lci year if model_n ==  6 & type_m ==  "OLS", 
            lcolor(black%80) )	
			,
		legend(order(1 "Benchmark PMT "		 			
                    3 "Updated PMT – employment status"					
                    5 "Dynamic " ) pos(6) row(2))
        ytitle("Exclusion Error")
        xtitle("Year")
        xlabel(2019(1)2021)
        ylabel(15(5)55)
        name(A, replace)
    ;
    #delimit cr
    graph export "FA4_performance_a_gainloss_`m_type'_bst.png", as(png) replace width(1200)
    graph export "FA4_performance_a_gainloss_`m_type'_bst.pdf", as(pdf) replace 
    graph export "FA4_performance_a_gainloss_`m_type'_bst.svg", as(svg) replace

    #delimit ;
    tw (connected inclusion_err year if type_m ==  "Baseline", 
            mcolor(red%80) msym(S) lcolor(red%80))
			(rcap inclusion_err_uci inclusion_err_lci  year if type_m ==  "Baseline", 
            lcolor(red%80))
        (connected inclusion_err year if type_m ==  "labour", 
            mcolor(ltblue) msym(0) lpattern(-) lcolor(ltblue))
			(rcap inclusion_err_uci inclusion_err_lci year if type_m ==  "labour", 
            lcolor(ltblue) lpattern(-) )
        (connected inclusion_err year if model_n ==  6 & type_m ==  "OLS" , 
            mcolor(black%80) msym(D) lcolor(black%80)) 
			 (rcap inclusion_err_uci inclusion_err_lci year if model_n ==  6 & type_m ==  "OLS" , 
            mcolor(black%80) msym(D) lcolor(black%80)) 
			,
        legend(order(1 "Benchmark PMT "		 			
                    3 "Updated PMT – employment status"					
                    5 "Dynamic " ) pos(6) row(2))
        ytitle("Inclusion Error")
        xtitle("Year")
        xlabel(2019(1)2021)
        ylabel(15(5)55)
        name(B, replace)
    ;
	
    #delimit cr
    graph export "FA4_performance_b_gainloss_`m_type'_bst.png", as(png) replace width(1200)
    graph export "FA4_performance_b_gainloss_`m_type'_bst.pdf", as(pdf) replace 
    graph export "FA4_performance_b_gainloss_`m_type'_bst.svg", as(svg) replace



    #delimit ;
    tw (connected crra year if type_m ==  "Baseline", 
            mcolor(red%80) msym(S) lcolor(red%80))
			(rcap crra_uci crra_lci year if type_m ==  "Baseline", 
            lcolor(red%80))
        (connected crra year if type_m ==  "labour", 
            mcolor(ltblue) msym(0) lpattern(-) lcolor(ltblue))
		        (rcap crra_uci crra_lci year if type_m ==  "labour", 
             lcolor(ltblue))	
        (connected crra year if model_n ==  6 & type_m ==  "OLS" , 
            mcolor(black%80) msym(D) lcolor(black%80)) 
		        (rcap crra_uci crra_lci year if model_n ==  6 & type_m ==  "OLS" , 
            lcolor(black%80))	
			,
        legend(order(1 "Benchmark PMT "		 			
                    3 "Updated PMT – employment status" 					
                    5 "Dynamic " ) pos(6) row(2))
        ytitle("Social Welfare (CRRA: {&rho} = 3)")
        xtitle("Year")
        xlabel(2019(1)2021)
        name(C, replace)
    ;
    #delimit cr
    graph export "FA4_performance_c_gainloss_`m_type'_bst.png", as(png) replace width(1200)
    graph export "FA4_performance_c_gainloss_`m_type'_bst.pdf", as(pdf) replace 
    graph export "FA4_performance_c_gainloss_`m_type'_bst.svg", as(svg) replace

    #delimit ;
    tw (connected hh_benefits year if type_m ==  "Baseline", 
            mcolor(red%80) msym(S) lcolor(red%80))
		(rcap hh_benefits_uci hh_benefits_lci year if type_m ==  "Baseline", 
             lcolor(red%80))	
        (connected hh_benefits year if type_m ==  "labour", 
            mcolor(ltblue) msym(0) lpattern(-) lcolor(ltblue))
		        (rcap hh_benefits_uci hh_benefits_lci year if type_m ==  "labour", 
            lcolor(ltblue))	
        (connected hh_benefits year if model_n ==  6 & type_m ==  "OLS" , 
            mcolor(black%80) msym(D) lcolor(black%80))
		        (rcap hh_benefits_uci hh_benefits_lci year if model_n ==  6 & type_m ==  "OLS" , 
            lcolor(black%80))	
			,
        legend(order(1 "Benchmark PMT"		 			
                    3 "Updated PMT – employment status" 					
                    5 "Dynamic" ) pos(6) row(2))
        ytitle("Per-household Transfer ($ USD PPP)")
        xtitle("Year")
        xlabel(2019(1)2021)
        name(D, replace)
    ;
    #delimit cr
    graph export "FA4_performance_d_gainloss_`m_type'_bst.png", as(png) replace width(1200)
    graph export "FA4_performance_d_gainloss_`m_type'_bst.pdf", as(pdf) replace 
    graph export "FA4_performance_d_gainloss_`m_type'_bst.svg", as(svg) replace

 }

